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^ . ABSTRACT 

' Aims. In this work we study the influence of the UV radiation field of a massive star on the evolution of a disklike mass of gas and 

pL^ , dust around a nearby star This system has similarities with the proplyds seen in Orion. 

Methods. We study disks with diff'erent inclinations and distances from the source, performing fully 3D numerical simulations. We 
use the YGUAZU-A adaptative grid code modified to account for EUV/FUV fluxes and non-spherical mass distributions. We treat 
H and C photoionization in order to reproduce the ionization fronts and photodissociation regions observed in proplyds. We also 
incorporate a wind from the ionizing source, in order to investigate the formation of the bow shock observed in several proplyds. We 

■ examine density and Ha maps, as well as the mass loss rates in the photoevaporated winds. 

' Results. Our results show that a photoevaporated wind propagates from the disk surface and becomes ionized after an ionization 

\ front (IF) seen as a bright peak in the Ha maps. We follow the development of an HI region inside the photoevaporated wind which 

corresponds to a photodissociated region (PDR) for most of our models, except those without a FUV flux. For disks that are at a 
distance from the source d >0A pc, the PDR is thick and the IF is detached from the disk surface. In contrast, for disks that are closer 
to the source, the PDR is thin and not resolved in our simulations. The IF then coincides with the first grid points of the disk that are 
Q I facing the ionizing photon source. In both cases, the photoevaporated wind shocks (after the IF) with the wind that comes from the 

?-H , ionizing source, and this interaction region is bright in Ha. 

jj^ ■ Conclusions. Our 3D models produce two emission features: a hemispherically shaped structure (associated with the IF) and a 

■ detached bow shock where both winds collide. A photodissociated region develops in all of the models exposed to the FUV flux. More 
importantly, disks with diff'erent inclinations with respect to the ionizating source have relatively similar photodissociation regions. 
If the disk axis is not aligned with the direction of the ionizing photon flux, the IF displays moderate side-to-side asymmetries, in 
qualitative agreement with images of proplyds, which also show such asymmetries. The mass loss rates are ~ 10~^ M© yr~^ for face-on 

■ disks, and 5 x 10~^ M© yr~^ for inclined disks at distances from 0.1 to 0.2 pc from the ionizing photon sources. 

\^ ' Key words. H II regions - hydrodynamics - stars: winds, outflows - circumstellar matter 

o ■ 
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, 1 . Introduction different parts of the object. Proplyds present a cometary shaped 

■ feature, very strong in Ha. Several of them show another bow 
^ The ultraviolet radiation from a massive star in an environment ^^aped feature, detached from the cometary structure itself, vis- 
^ of low mass star formation produces several interesting phenom- .^^^ 5qq7 ^ j^^^^^^ j ^^^^^^ ^^^^ 

^ ena, among them the appearance of the so-called proplyds (an ^j^^ ^^^^ ^^^^^ ^j^^ ionization front (see Figure D. In 
. ^ acronym that stands for PROto-PLanetarY Disk; see O Dell et 
^ , al. 1993 and O'Dell and Wen 1994). Proplyds are young stel- 



Oh: 



some objects, the disks are seen in emission in [01] 6300 A. 

^ 1 . /^zc^r^\ VI . 1 J • r . .1 . It is also possible to see that the optical lines are double or even 

, lar obiects (YSO) with a cometary- shaped ionization front that iji- i -.i .i^r.i ij 

rri ^ . ^ J J u triple peaked. For example, m the central parts or the proplyd 

^ ■ contain a central, low mass young star surrounded by an ac- . of-^ /t ^^^x tt i fi . i i 

■"■ ^-i-i 1^ Ij -j 167-317 (LV2), Hof has three components: a low velocity com- 

cretion disk, a photoevaporated wind, and, m some cases, a . i • ' i- i i v r o^m -i i • i 

1 il - 1 .r, rx., 1 J • .1 • u ponent, peaking at a radial velocity or ~ 30 km s \ and two high 
monopolar or a bipolar outflow. The proplyds m the Orion neb- i . . a . ^f^f^^ -i j nc ^ 

1 /AT/io\ a ^ ^ A 1 J ^- 1 • velocity components, centered at ~ 100 km s 'and 75 km 

ula (M42) were first report ed as unresolved optical emission /j y^g^Qj^ J^^Qg ^\ I | 2Q05 b 
fe atures by Laques & Vidall ([l979i) and then as radio sources ^ ^' 

bylChurchwefl et al. Cll987h. Using the Hubble Space Telesc ope Johnstone et alj ([l998h explained the observations with a 
(HST), O'Dell & We^ 119941) and lO'Defl & Wong (|199^ re- model in which the radiative field and the wind from 6^ Ori C in- 
ported the discovery of tens of proplyds around the Trapezium teract with a low mass YSO, generating a photoevaporated wind, 
cluster. an ionization front (IF) and a bow shock. They also proposed 
The proplyds present strong optical e mission lines such the existence of a photodissociated region (PDR), in which the 
as Hof, [OIII] 5007 A, [ 01] 6300 A (Bally et all, |200J ) photoevaporated wind is neutral (the H is neutral). They sepa- 
and also sorne ultra violet (iHenney et al.l |2002) and IR lines rate their models in two classes that depend on the amount of 
(iTakami et al.L l2002h . With HST high spatial resolution images EUV and FUV flux arriving at the proplyd's disk surface. The 
and spectra, it was possible to see that the emission comes from EUV model, valid for a YSO near to or very far from 6^ Ori C 
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Fig. 1. A cartoon of a star-disk system being photoevaporated by 
the UV radiation field from an O type star. The major features are 
highlighted: the proplyd head and tail which traces the ionization 
front from the direct ionization field coming from the 0-star and 
from the diff'use radiative field, respectively, and the interaction 
region between the photoevaporated flow and the wind from the 
0-star. 



(d < 10^^ cm or J > 10^^ cm), presents a very thin FDR. In this 
case, the EUV flux determines the photoevaporated mass loss 
rate. On the other hand, for the FUV model, the PDR is thick 
and can contain a shock fr ont with in it. Here, the FUV flux de- 
termines the mass loss r ate.lStorzer & HollenbachI (1 1 999l) further 
improved the models by lJohnstone et al.. (.1998.) including results 
from PDR codes. 

Some axisymmetric numerical simulations were also 
done in order to reproduce the characteristics of proplyds. 
iGarcia-Arredondo et al.l (12001 1) carried out axisymmetric simu- 
lations of the interaction of the stellar wind from 6^ Ori C with 
the photoevaporated proplyd flow. They took available observa- 
tional data to constrain stellar parameters such as the wind den- 
sity and velocity and the ionization front parameters. Their nu- 
merical calculations also assume a hemispherical proplyd head 
and a cylindrical tail. Studying the interaction between these two 
winds, they have successfully reproduced the arc emission for 
the proplyds near O^^C. 

'Richli ng & Yorkd (l2000h performed 2D, axisymmetric hy- 
drodynamical simulations of photoevaporating disks including 
both ionizing (hv > 13.6 eV) and dissociating (6 eV < hv < 
13.6 eV) radiation. In their models, disk structures are formed 
through collapse simulations o f 1 and 2 M© molecular clumps 
(lYorke & Bodenheimeii Il999l) . which are then exposed to the 
radiation field by switching on the external UV radiation field 
in the calculation. They studied the eff'ects of distance from the 
UV photon source on emission line maps and the eff'ect of the 
presence of a spherical wind from the proplyd star, which pro- 
motes the appearance of collimated, bipolar microjets. Instead 
of calculating the dissociation of the H2 molecule, they followed 
the ionization of C. They took into account photoelectric heating 
and cooling by fine-structure lines such as [CII] 158 yum, [01] 
63 yum and [O I] 145 jim. An important point in their work is the 
inclusion of the diff'use radiation field, which is responsible for 
the tails of the proplyds (see Figure [TJ. 

In this work we present the first fully three-dimensional nu- 
merical simulations of disks exposed to FUV and EUV radiation 




Fig. 2. Schematic view of the initial setup of the simulations. The 
cartoon shows also the distance from the disk to the ionization 
front, rjF, and the definition of the width of th e ioniz ed structure, 
r^. We follow the notation of I Johnstone et al.l (1 1998b (see Figure 
7 in that paper). 

Table 1. Computed models 



Model 


Spuv 


d 






T^ 




(Xl049 s-1) 


(pc) 




(km s-^ 


(K) 


Mia 





0.1 


0° 





1000 


Mlb 





0.1 


0° 


100 


1000 


M2a 


1.78 


0.2 


0° 


100 


1000 


M2b 


1.78 


0.1 


0° 


100 


1000 


M2c 


1.78 


0.02 


0° 


100 


1000 


M2d 


1.78 


0.1 


0° 


100 


3 000 


M3a 


1.78 


0.2 


45° 


100 


1000 


M3b 


1.78 


0.1 


45° 


100 


1000 


M4a 


1.78 


0.2 


75° 


100 


1000 


M4b 


1.78 


0.1 


75° 


100 


1000 



^ Inclination of the disk, in degrees, with respect to the z-axis. 

^ Velocity of the wind from the O star. 

^ The temperature associated with the PDR. 



fields. As discussed in I Johnstone et al.l (Il998l) , the diff'use ioniz- 
ing and dissociating radiation field is important for determining 
the shape of the ionization front behind the (directly) illuminated 
disk surface. Our present simulations do not include the diff'use 
field, resulting in the tails that do not have the correct morphol- 
ogy (see Cerqueira et al. 2006b). We concentrate on obtaining a 
description of the head of the proplyd flow, and study the efl'ect 
of difl'erent orientations of the disks with respect to the imping- 
ing UV photon field (and the stellar wind). 

The present paper is organized as follows. In ^ we de- 
scribe the code and the simulations. In ^ we present our results. 
Finally, in ^ we draw our main conclusions. 

2. The numerical method, the numerical setup, and 
the simulated models 

The 3D numerical simulations ha ve been carried ou t with the 
YGUAZU-A adaptive grid code (iRaga et al.L l2000l 12002) us- 
ing a 5 -level binary adaptive grid. The YGUAZU-A code in- 
tegrates the gasdynamic equations employing the flux vector 
splitting scheme of van Leer (1982) together with a system 
of rate equations for atomic/ionic species. In our simulations, 
we consider 4 species: HI, HII, CI and CII. This code has 
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Fig. 3. Density xz-midplane stratifications (left) and Ha maps (right) obtained from models Mia (top) and Mlb (bottom) which has 
only an EUV radiation field for r = 2 000 yr. The scales (given by the top-bars) are in g cm"^ for the density (left), and normalized 
to one for the Ha maps (right). The coordinate axes are in AU. 



been extensively employed for simulating diff^e rent astrophys- 
ical flows such as jets (Masciadri et al., 2002; C erqueira et all 
l200 6a), interacting winds dGo nzalez et al., 2004), photoevapo- 
rating clumps ( Cerqueira et al.t. 2006bi) and supernova remnants 
(IVelazquez eTaP. l2001aL l2004l) . It was also tested with laser 
generated p lasma laboratory experiments (ISobral et al.i 120001: 
iRaga etal.L f2001: Velazquez et al, 2001^. 

In this work, instead of solving an energy equation, we pre- 
scribe a temperature law, given by 



T = (Ti- T2) • xhii + T2 ' xcu + • (1 - xcii), 



(1) 



where Ti = 10 000 K, is the characteristic temperature of an H II 
region, T2 = 1 000 K, the typical PDR temperature, ^3 = 10 K, 
the temperature of the molecular gas, xuii is the hydrogen ioniza- 
tion fraction and xcii is the carbon ionization fraction (xcii = 1 
when all the CI is CII). This prescription is justified if the ther- 
mal equilibrium time scale is much sm aller than the dynamical 
time scale (Leflo ch & Lazarelfl Il994l) . which is the case here, 
for both the ionized and the molecular gas. It is clear from equa- 
tion ([T]), that we have possible temperatures ranging from 10 K 
to 10^ K. In order to study the dependence of the PDR geome- 
try with temperature, we also compute models in which we set 
7^2 = 3 000 K. 

Following iRichling & York3 (l2000b we do not treat the dis- 
sociation of H2, but consider that the dissociation front and the 
carbon ionization front coincide. We then solve rate equations 
for hydrogen and carbon taking into account photoionization, 
collisional ionization, radiative and dielectronic recombination. 

We model the O star as a point source located outside the 
computational domain. This source emits EUV photons (with 



hv > 13.6 eV) at a rate ^euv = 7.2 x 10^^ s"^ and FUV photons 
(with 6 eV < hv < 13.6 eV) at a rate 5fuv = 1-78 x 10^^ s'^ 
(or zero, in order to study the efl'ects of the EUV radiation field 
only; see Obelow). 

For the transfer of FUV photons, we consider the optical 
depth due to the ionization edge of the C I photoionization cross 
section, and to the absorption of photons by dust: 



tfuv = 1.22 X IQ-^^cm^Nci + 1.15 x IQ-^^cm^ Nh 



(2) 



where Nci is the C I column density integrated along the path 
of the ray. The second term (representing the FUV dust absorp- 
tion) is calculated assuming that the optical extinction is given 
by Ay = Nhi/2 x 10^^ cm"^ (appropriate for a standard grain 
distribution). We then multiply this number by 2.4 to obtain the 
FUV extinction. This is consistent with empirically derived op- 
tical/UV extinction curves (see, e. g., Cardelli et al. 1988). We 
are also implicitly assuming that there is no dust in the region 
where H is fully ionized. 

Our simulations include a flat, neutral disk, with a constant 
thickness H = 15.6 AU and an outer radius of 50 AU. The disk is 
positioned at distances of 0.02, 0.1 and 0.2 pc from the ionizing 
photon source. With this particular choice of parameters, we can 
in principle ensure that the photoevaporated flow of the near- 
est disk would be EUV-dominated, while for the more distant 
one, FUV-dominated. Using previous analytical calculations of 
I Johnstone et al.l (Il998l) and the disk radius given above, we can 
roughly estimate the minimum distance from the ionizing source 
in order to have an FUV-dominated flow as dmin ~ 0.1 pc. The 
disk radius that would support an FUV-dominated flow is, for 
d = 0.1 pc, ^ 30 AU, and ^ 120 AU for d = 0.2 pc 
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Fig. 4. Density xz-midplane stratifications for model Mlb, which 
has only an EUV radiation field, from ^ = 80 yr (top) to t = 120 
yr (bottom), showing the interaction between the stellar wind 
with the photoevaporated flow. The density scale (given by the 
top-bar) is in g cm"^. The coordinate axes are in AU. 



(see equations 15 and 16 in Johnstone et al 1998). We choose to 
simulate disks at distances from the ionizing photon source that 
are i) at a distance that is suitable for an EUV-dominated flow 
(d = 0.02 pc), ii) close to the limiting cases for which we should 
expect an FUV dominated flow to occur (d = OApc^r^ = 50AU), 
and iii) FUV-dominated flow (d = 0.2pc). In this last case, the 
disk radius = 50 AU is found to obey the < 122 AU con- 
dition, and the system would develop an FUV-dominated flow. 
The only caveat is that, as pointed out by Johnstone et al. (1998), 
there are several sources of uncertainties in evaluating these (an- 
alytical) numbers. 



The density of the disk is assumed to be constant in both the 
vertical and radial directions. We start the simulations assuming 
pressure equilibrium between the PDR region and the disk itself: 



nd = no(T2/T3), 



(3) 



where T2 and ^3 are the PDR and disk temperatures (see 
Equation [T]) and no is the number density at the base of the pho- 
toevaporated flow, which can be estimated as: 



no - No/rd. 



(4) 



for a FUV dominated flow. Here, A^^ is the column density of the 
PDR attained for tfuv ~ 1 - 3, that is of the order of 10^^ cm"^ 
(iJohnstone et al.Lll998h . Using rd = 50 AU, we obtain no ^ 10^ 
cm~^. We therefore set nd = 10^ cm"^. For an EUV dominated 
flow, the density is found by equating the EUV flux to the re- 
combinations in the ionized flow: 



S EUV ^ 2 



(5) 



where ar = 2.6 x 10"^^ cm^ s"^ is the recombination coefficient 
for Hydrogen at lO^K,d is the distance from the disk to the ion- 
izing photon source, and S euv = 7.2 x 10^^ s"^ and = 50 AU 
in all of the models presented in this paper. We should note that 
the factor of 3 on the right hand side of equation (|5]) is a result of 
assuming a spherical divergence for the photoevaporated wind. 

We then calculate the densities no for diff'erent distances and 
find no = (0.98 - 5.0) x 10^ cm"^ for d = (OA - 0.02) pc. For 
all cases, we treat the disk as an infinite reservoir, since at each 
time step we reset its initial conditions, as described above. 

We assume the presence of a O.2M0 star at the disk center 
and consider the gravitational force that it exerts on the flow (by 
including the appropriate source terms in the momentum and 
energy equations). A Keplerian rotation velocity law (consistent 
with the 0.2 M© central star) is imposed on the disk material. 
The gravitational force is important as it inhibits the formation 
of a photoevaporated wind in the central regions of the disk. We 
expect a photoevaporated flow to start beyond the gravitational 
radius, which can be estimated as = GM^Ia^, where a is 
the local sound speed. For FUV dominated flows, a ^ 3 km 
s"\ so that Vg = 2.9 x 10^"^ cm, or ^ 20 AU. The inhibition of 
the photoevaporated wind in the central regions of the disk can 
be seen in our FUV dominated flow simulations (which have 
rg ^ 20 AU), but not in our EUV dominated flow simulations 
(which have rg ^ 2 AU, smaller than our numerical resolution). 
The size of the remnant disk, which is known to be typically 
smaller than the gravitational radius, will not be investigated in 
this work, since we are treating the disk as an infinite reservoir 
(IJohnstone etal.Lll998h . 

For some models, we add a plane-parallel wind (flowing par- 
allel to the ionizing/dissociating external radiation field, with ve- 
locity Vw = 100 km s"^ and with number density n^ = 500 
cm"^), which represents the wind from the 0-type star. The wind 
velocity of 0^ Ori C is estimated to be of the order of 500 to 
1 650 km/s (iBallv et al .. 1998), or even larger (up to about 3 600 
km s"\ Walbom & Nichols 1994), which is by far larger than 
the adopted velocity here. On the other hand, the ram pressure 
estimated for both the stellar wind and the photoevaporated flow 
for proplyds in the Trapezium region is ~ 3 x 10"^ dyn cm"^ to 
5x10"^ dyn cm"^ for the proplyds closer to 0^ Ori C (at distances 
^ 0.01/7C; see, e.g., Garcia- Arredondo et al 2001). Therefore, 
our particular choice of parameters for the stellar wind give us 
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Fig. 5. Distances from the disk to the interaction region as a function of time are depicted for model Mlb, which has only an 
EUV radiation field: <id-iR,c (full line) are the calculated distances, where stellar and the photoevaporated wind ram pressures match; 
<^d-iR,m (dotted-line) are the distances as measured from the simulations (along the x-axis, with y = z = 4.9x 10^^ cm). The distances 
are in units of 10^^ cm (left-axis). Also shown is the disk mass loss M (dashed-line), in units of 10"^ M© yr"^ (right axis). The M 
has been calculated by computing the quantity pA • v on the faces of a cube centered at the disk center, with a = 3rd. 



a ram pressure comparable to that estimated for a group of pro- 
plyds in the Trapezium cluster. In particular, with the adopted 
parameters, we have Psw = 8.3 x 10"^ dyn cm"^. 

The initial setup of the simulations is shown in Fig. [2l The 
flat, dense structure with which we model the circumstellar disk 
is placed on the }7z-plane, at x = 1 203 AU. The disk axis, which 
we define here as the axis perpendicular to the disk plane, and 
that contains the disk center at (x,y,z)= (1203,334,334) AU, co- 
incides with the X-axis for an inclination angle ^ = 0° . For larger 
inclination angles, the disk axis rotates clockwise, with respect 
to the X-axis, as depicted in Figure [2l The remaining volume of 
the computational domain is filled with an ionized environment 
of uniform density Ua = 500 cm"^. All models except Mia have 
an initial wind of v^ = 100 km s"^ and density n^ = 500 cm"^, 
entering the computational domain through the x = bound- 
ary. The ionizing/dissociating photon source is placed along the 
X-axis, outside the computational domain to the left. 

In Table [T] we present all the relevant parameters used in our 
models: the dissociating photon rate, S fuv (in units of 10"^^ s"^; 
second column), the distance from the disk to the external star 
(in pc; third column), the inclination angle 6 between the disk 
and the x axis (fourth column), the velocity of the wind from 
the external star (in km s"^; fifth column) and the temperature 
of the PDR region (in K; sixth column). All the models have an 
ionizing photon flux S euv = 7.2 x 10"^^ s"^ (similar to the values 
used by Richling & Yorke 2000) and = 50 AU. 

The simulations have been computed in a 5 -level, binary 
adaptive grid with a maximum resolution (along the three axes) 
of 5.2 AU for all models. The computational domain has a size 

(Xmax,ymax,Zmax) =(1 337,668,668) AU. 



3. Results and Discussion 

In Fig. O we show the xz-midplane density stratifications (left) 
and Hof maps (right) for model Mia (top) and model Mlb (bot- 
tom) at t = 2 000 years. The models present only an EUV ra- 
diation field from the ionizing source (at a 0.1 pc distance from 
the disk). The diff'erence between the two models is the presence 
of a stellar wind in model Mlb, as described in Table [T] The 
disk axis has no inclination with respect to the direction of the 
impinging photons (i. e., ^ = 0° see Figure |2] and also Table [TJ. 

The Ha emission in model Mia, is restricted to a region 
near the disk surface (and facing the ionizing source), and the 
edge of the shadow on the non-illuminated side of the disk. For 
model Mlb, we have Ha emission close to the disk (with mor- 
phology and intensity similar to those seen in model Mia), and 
also fainter emission regions distributed along the interaction re- 
gion between the stellar wind and the photoevaporated flow. It 
is known that some proplyds in Orion show faint Ha arcs fac- 
ing 6^ Ori C (Bally et al., 1998), which correspond to the stel- 
lar wind/photoevaporated flow interaction emission found in our 
models. 

In Figure |4] we show xz-midplane density maps for model 
Mlb at t = 60, 80, 100 and 120 yr, in order to illustrate the ini- 
tial evolution of such an interaction. At ^ = 60 yr, the stellar wind 
and the photoevaporated flow start to interact at x ^ 668 AU. At 
later times, the interaction region moves towards the disk, and 
stabilizes around the position of the point of ram pressure bal- 
ance between the stellar wind and the photoevaporated flow. We 
can also see this efl'ect in Figure [51 where we show the photo- 
evaporated wind mass loss rate (M; dashed line), the distance 
from the disk to the interaction region i) as measured directly 
from the xz-midplane density maps (<id-iR,m; dotted-line) and ii) 
as calculated (Jd-iR,c; full-line) by balancing the ram pressures 
from the stellar wind (Psw) and from the photoevaporated flow 
(Ppef)- The expression for <id-iR,c, that takes into account the log- 
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Fig. 6. Density xz-midplane stratifications (left) and Ha maps (right) for models M2a (top), M2b (middle) and M2c (bottom) at 
t = 2 000 yr. The scales (given by the top bars) are in g cm"^ for the density (left), and normalized to one for the Ha maps (right). 
The coordinate axes are in AU. 



arithmic outward acceleration expected for the photoevaporated 
wind is given by (see Henney et al. 1996): 



^^d-IR,c - ^0 ■ 



l+41n|^ 



1/4 



where is the disk radius, and Rq: 



j Mail 



(6) 



(7) 



In this equation, an is the sound speed of the photoionized flow, 
and Psw is the ram pressure of the steflar wind, which is Psw = 
8.3x10"^ dyn cm"^ for our chosen parameters (see ^2]). In Figure 
\5\ we see that both the measured and the calculated distances are 
comparable. 

In Fig. O we show the xz-midplane density stratifications 
(left) and Ha maps (right) for models M2a (top), M2b (mid- 
dle) and M2c (bottom), ait = 2 000 years. The models have both 
EUV and FUV radiation fields (from the stellar source). The dis- 
tance from the disk to the ionizing source is 0.2 pc, 0.1 pc and 



0.02 pc for models M2a, M2b and M2c, respectively; see Table 
[TJ The density of the disk is determined from equations ^ and 
Q for the FUV dominated models M2a and M2b, and equations 
0]) and ([5]) for the EUV dominated model M2c. The disk axis is 
aligned with the x-direction (6 = 0° see Figure [2] and also Table 
[T]), and the stellar wind is present in the three models. 



We see that a high density region shrouds the disk in mod- 
els M2a and M2b. This region is related to the photodissociation 
region exp ected from the analytic al models of Jo hnstone et al.l 
(Il998h and IStorzer & HoUenbachl fl999). It has a temperature 
72 - 1 000 K (see equation [T] and Table [T]). This region is com- 
posed of C II and H I. We also note that the environment (the left 
part of the computational domain) and the photoevaporated wind 
outside the PDR have a temperature of about 10"^ K, indicating 
that they are photoionized. The photodissociation region sepa- 
rates the disk from the ionized photoevaporated wind. In model 
M2c (in which the disk is at 0.02 pc from the ionizing photon 
source) the photodissociation region is limited to a thin shell in 
contact with the surface of the disk. 
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Fig. 7. Density xz-midplane stratifications superimposed on velocity field vectors (left panels) and an xz-map of the Mach number 
(right panels), M = v/a (a is the local sound speed) for models M2a (top) and M2b (bottom), at ^ = 200 yr. The snapshots are a 
closer view of the star-disk system being photoevaporated. The distance (in AU) to the origin of the coordinate system is depicted 
for the X and z axes in the left panels. The right-hand side bars on each panel indicate the density scale (in g cm"^; left panels) and 
the Mach number (right panels). The labels in the top-rigth panel indicates the region where we have a supersonic outflow (1 and 
3), a subsonic outflow (2), and a subsonic inflow (4). We also note that region (1) is neutral, at a temperature T2 = 10^ K, whereas 
(3) is ionized, at a temperature Ti = 10^ K. 



The ionization front radius (as defined in Figure O is found 
to grow from rjp = 2.4rj = 120AU at t = 200 yr to rjf = 
3.8rj = 190AU at ^ = 2 000 yr for model M2a, while for model 
M2b, rjF = 1.6rj = 80 AU (constant throughout the simulation). 
In models M2a and M2b, the outer radius of the photodissocia- 
tion region (rw, defined in Figure 2) grows with time. This lateral 
expansion occurs at a velocity of ^ 0.1 km s"\ at ~ 1/10 of the 
~ 1 km s"^ outwards flow velocity of the gas within the pho- 
todissociation region. This indicates that the ionization front in 
these models is pressure confined by the base of the photoevapo- 
rated flow. The ratio between the cross section radius and the 
axial standofl' distance rjF (see Figure 2) has values r^/riF = 1.4 
and 1.9 for m odels M2a and M2b (at t = 2 000 yr), respectively. 
Interestingly. I Johnstone et al] (1 19981) predict = l.SrjF from 
their simplified, analytic model. Models M2a and M2b have 
detached, almost hemispherical ionization fronts. On the other 
hand, in model M2c the ionization front lies very close to the 



surface of the disk (with just one pixel of separation between 
them ). It is interesting to no te that observations of proplyd 170- 
337 (iJohnstone et al.L Il998h . that is situated at 0.1 pc from 6^ 
Ori, show that its IF radi us is approxima tely 80 AU. From the 
analytical calculations of Johnstone et al. (1998) for FUV dom- 
inated flows, riFi4 - OA(e^/ frS 4g)d^^, where riFi4 is the radius 
of the IF in units of 10^^ cm, Jn is the distance from the source 
in units of 10^^ cm, ^"49 is the ionizing photon rate in units of 
10"^^ s"\ /r is the fraction of FUV photons absorbed by recom- 
bination in the ionized region of the flow , and e i s a dim ension- 
less parameter o f order unitv (see Johnstone et al.l[l998l) . Taking 
(e^/frS 49) = 16 (IJohnstone et al.Lll998l) . and di7 = 3.08, we ob- 
tain riF ^100 AU, which is consistent with the values obtained 
from models M2a and M2b. 

In Figure E] we show the velocity field in a region close to the 
accretion disk superimposed on density (left panels) and Mach 
number maps (right panels), M = v/a, where a is the local sound 
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speed and v is the local velocity of the gas, for models M2a (top 
panels) and M2b (bottom panels), ait = 200 yr. From Figure H 
we see that the photoevaporated flow has four difl'erent dynam- 
ical regimes, indicated in the top-right panel of Figure |7]by the 
labels 1 to 4, which refer to: 1) a supersonic neutral outflow, 2) 
a subsonic, neutral outflow, both of them inside the PDR region 
and separated by a shock; 3) a supersonic, ionized outflow out- 
side the PDR region and 4) a subsonic, neutral inflow, localized 
around the disk center (at r < r^). 

As we have stated before in ^ the disk radius adopted in 
our models is larger than the gravitational radius. For instance, 
for the parameters of the M2 models (M2a, M2b and M2c), we 
have rg ^ 20 AU. We can see from the velocity field maps that 
between disk radial distances r = (the disk axis) and r ^ ±3Az 
(where Az = 5.2 AU is the spacing of the highest resolution 
grid, see ^2]), we do not have an outflow. For r > 20 AU, the 
outflow starts supersonically from the disk surface, with Mach 
numbers 1 . 1 (region 1 in Figure [7]), and reach an internal shock 
(before the IF; see Figure|7]) with Mach numbers 1.8 (note that 
the PDR has a sound velocity aj = 3.2 km s"^). Between this 
internal shock and the outer boundary of the PDR, the outflow 
is subsonic, with Mach numbers ^ 0.4-0.6 (region 2 in Figure 
[7]). Outside the PDR and before reaching the interaction region 
between the two winds, we have a supersonic outflow with Mach 
numbers ^ 1.2-1.6 (with a sound velocity an = 14.2 km s"^). 

The structure obtained f rom models M2a and M2b is similar 
to the scenario proposed by I Johnstone et al.l (Il998b for the FUV 
dominated flow. From this we conclude that both the M2a and 
M2b models produce photoevaporated flows in the FUV domi- 
nated regime. 

In Figure [8] we show the xz-midplane density stratifications 
(left) and Ha maps (right) for models M3a (top) and M3b (bot- 
tom panels) at ^ = 2 000 yr. The models have an FUV -h FUV 
radiation field from the ionizing source, that is located at 0.2 pc 
(M3a) and 0.1 pc (M3b) distance from the disk. In both cases, 
the disks have a 45° of inclination with respect to the direction 
of the impinging photons (^ = 45° see Figure |2] and Table [T]). In 
Figure [9] we present the same results, but for models M4a and 
M4b (with = 75°). As can be seen in Table [H models M2a, 
M3a and M4a diff'er only by the values of the disk inclination. 
The same is true for the set of models M2b, M3b and M4b. 

In the inclined disk models (models M3a, M3b, M4a and 
M4b, see Figures 8-9) the internal structure of the PDR region is 
the same as in the M2a and M2b models (with ^ = 0°, see Figure 
[7]and Table 1). There is a neutral supersonic outflow (with M ^ 2 
in all models) near the illuminated surface of the disk. For inner 
disk radial distances (from the disk axis), there is a gas inflow. 
In Figure [TOl we show velocity vectors superimposed on density 
maps for models M3a (top-left), M3b (top-right), M4a (bottom- 
left) and M4b (bottom-right). As in the M2a and M2b models, 
we again have an internal shock that separates the supersonic 
from the subsonic neutral flows, inside the PDR. Outside this 
internal shock we have the outer edge of the PDR region, that 
separates the subsonic, neutral flow from the external, ionized 
supersonic flow. These structures are characteristic of an FUV 
dominated flow, as discussed previously. 

The Ha emission in the region close to the ionization front of 
models M3 and M4 (Figures [8] and O has a general morphologi- 
cal resemblance to the Ha morphology of the "face-on" models 
(M2, see Figure [6]). However, the Ha maps of the inclined disk 
models show clear side-to- side asymmetries, which are more 
pronounced for the models with ^ = 45° (M3a, b, see Table 1 
and Figure (8]) than for the ones with = 75° (M4a, b, see Figure 
i). 



From the t = 2 000 yr time-frames of models M3 and M4 
(see Figures 8 and 9), we have calculated the axial stand-ofl' dis- 
tance (rjF) and width of the cross section presented by the ion- 
ization front (r^, see Figure 2). We obtain the following results: 

- M3a: rjF = 3.6rj = 180 AU; = 4.4rj = 220 AU, 

- M3b: riF = 2.0rj = 100 AU; = 2.4rj = 120 AU, 

- M4a: rjF = 2.8rj = 140 AU; = S.lr^ = 185 AU, 

- M4b: riF = 1.9r^ = 95 AU; = 2.1r^ = 105 AU, 

with r^/rjF < 1.3 for all models. These results illustrate the re- 
duction in the size of the photodissociation region that is ob- 
tained for increasing values of the disk inclination angle 6. 

In Figure [TT] we present the mass loss rate as a function of 
time, for all the simulated models except model M2d. The mass 
loss rate for models Mia and Mlb (without an FUV radiation 
field; see Table [B shows a small time variability (as mentioned 
above), but is almost independent of the presence of the stellar 
wind. Also, the mass loss rate has approximately the same value 
even if the distance from the disk to the ionizing photon source 
changes by a factor of two, even for the are FUV dominated 
models, as is the case for those with a given inclination angle 
(with respect to the ionizing photon source) as M2a and M2b 
(0 = 0°), M3a and M3b (0 = 45°), and M4a and M4b (0 = 75°). 
This is due to our simplified model, in which the PDR temper- 
ature is assumed to be position independent and with the same 
value for all models. Looking at a series of models with increas- 
ing values of the inclination angle (e. g., models M2a, M3a 
and M4a, which otherwise have identical paramenters, see Table 
1), we see that the mass loss rate as a monotonically decreasing 
function of 0. This result is to be expected, given the lower eff'ec- 
tive cross sections (to the impinging UEV/FUV radiation field) 
presented by the disk for increasing inclination angles. In partic- 
ular, if we take the mean value of M for models M2a and M2b 
(with e = 0), for models M3a and M3b (with 6 = 45°) and for 
models M4a and M4b (with = 75°), we find Mm2 = 1.4 x 10"^ 
Mo yr-\ Mm3 = 6.2 x 10"^ Mq yr-^MM4 = 4.8 x 10"^ M© yr-\ 
which implies Mm3 - 0AM mi and Mm4 - 03Mm2- 

In Fig. [121 we show the xz-midplane density stratifications 
(left) and Ha maps (right) for model M2d at t = 200 yr. This 
model has the same parameters as the M2b model (6 = 0°, d = 
0.1 pc; see Table© but a diff'erent PDR temperature T2 = 3 000 
K (see equation [T]). The ionization front radius is larger, by a 
factor of ^ 3, for the model M2d when compared with those 
with T2 = 1 000 K. The mass loss rate is 5 x 10"^ M© yr"^ for 
M2d (about 3.5Mm2). We have run equivalent simulations for 
models M3b and M4b, that is, a model with 6 = 45° and d = OA 
pc, but with T2 = 3 000 K, and a model with 6> = 75° and J = 0. 1 
pc, but with T2 = 3 000 K. For these models (not shown here), 
we found the same trend of larger mass loss rates and ionization 
front standoff' distances for larger values of T2. 

4. Conclusions 

This work describes the first fully three-dimensional numeri- 
cal gasdynamic simulations of the photoevaporation of neutral 
disks, subject to an interaction with the FUV/EUV photon ffux 
and the wind from an external O star. 

With this initial setup, we study the time-evolution of the 
photoevaporated ffow. We have calculated models in the FUV 
dominated regime, in which one obtains an ionization front (IF) 
close to the surface of the disk, and a photodissociation region 
(PDR) not resolved by our numerical scheme between the sur- 
face of the disk and the IF, and in the FUV dominated regime. 
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Fig. 8. Density xz-midplane stratifications (left) and Ha maps (right) obtained from models M3a (top) and M3b (bottom) for t = 
2 000 yr. The scales (given by the top bars) are in g cm"^ for the density (left), and normalized to one for the Ha maps (right). The 
coordinate axes are in AU. The distance between the ionized photon source and the disk is d = 0.2 pc and J = 0.1 pc for models 
M3a and M3b, respectively. In both cases, = 45° (see Figure 2). 



in which a PDR (well resolved in our numerical simulations) is 
formed between the surface of the disk and the IF. Both EUV 
and FUV dominated models develop an interaction region be- 
tween the photoevaporated wind and the wind from the hot star. 
In some cases, this interaction region lies within our computa- 
tional domain. For these cases, the interaction region lies close 
to the position of stellar wind/photoevaporated flow ram pressure 
balance. 

We have simulated models of disks with difl'erent inclina- 
tions with respect to the direction of the impinging EUV and 
FUV photon fluxes, and at difl'erent distances from the ionizing 
photon sources ranging from 0.02 pc to 0.2 pc. The model at a 
distance of 0.02 pc from the 0-star is clearly EUV dominated, 
since the IF is just above the disk surface. Models at distances 
from 0.1 pc to 0.2 pc show detached PDRs. These PDRs show 
structures, such as the presence of a neutral supersonic flow near 
the disk and, farther out but before reaching the IF, a subsonic 
neutral flow, that are consistent with the main features predicted 
by Johnstone et al. (1998) for a FUV-dominated flow. We note 
that these structures occur in our simulations even for the in- 
clined disk models. The typical size of the IF was found to be 
in agreen ient with analytical resu lts. In particular, the limits pre- 
dicted by I Johnstone et al.l (1 19981) for the ionization front radius, 
2.5rj < riF < 4rj, were recovered in most of the models, with 
or without disk inclination. 

We found that, for a given inclination of the disk, the 
mass loss rate seems to be independent of the distance to the 
FUV/EUV photon source (for the FUV-dominated flows). In our 
simulations this is due to the constancy of the PDR tempera- 
ture. For higher disk inclinations, we obtain smaller mass loss 



rates, which is a direct result of the smaller efl'ective cross sec- 
tion presented by the disk to the impinging FUV/EUV photon 
flux. Taking the mean values of the mass loss rate from our sim- 
ulations, we obtain Mm2 = 1.4x10"^ Mo yr"\MM3 = 6.2x10"^ 
Mo yr-iMM4 = 4.8 x 10"^ Mo yr"\ for models with 6 = 0°, 
^ = 45° and 6 = 15°, respectively. These values are surprisingly 
similar, in spite of the very difl'erent inclination angles used in 
the simulations. 

Our models predict Ha emission maps with moderate side- 
to-side asymmetries, which directly result from the misalign- 
ment of the axis of the accretion disk with respect to the direc- 
tion to the photoionizing source. Asymmetries are seen both in 
the emission associated with the IF and in the emission of the 
photoevaporated wind/expanding H II region interaction work- 
ing surface. It is clear that asymmetries in these two structures 
are indeed observed in many of the Orion proplyds (iBally et al.l 
2000). A detailed attempt to model a specific, asymmetrical pro- 
plyd with a 3D photoionizing disk simulation will be presented 
in a future paper. 
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Fig. 10. Density xz-midplane stratifications superimposed on velocity field vectors for models M3a (top-left), M3b (top-right), M4a 
(bottom-left) and M4b (bottom-right), at ^ = 2 000 yr. The distance to the origin of the coordinate system, in AU, is depicted for the 
X and z axes. The right-hand side bars indicate the density scale (in g cm"^). 
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Fig. 11. The mass loss rate as a function of time for models Mia and Mlb (top-left panel), M2a, M2b and M2c (top-right panel), 
M3a and M3b (bottom-left panel) and M4a and M4b (bottom-right panel). The caption on the top of each diagram indicates the 
model that a given curve type represents. The mass loss rate is in units of 10"^ M© yr~^ for models M2, M3 and M4 and in units of 
10"^ Mo yr"^ for models Ml. 




Fig. 12. Density xz-midplane stratification (left) and Ha map (right) for model M2d at ^ = 200 yr. The temperature of the PDR 
region in this model is 72 = 3 000 K. The scales (given by the top bars) are in g cm~^ for the density (left), and normalized to one 
for the Ha maps (right). The coordinate axes are in AU. 



